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Abstract 

In this research article, we demonstrate the application of the 
New Homotopy perturbation method to work out the 
nonlinear singular boundary value problems which often 
arise in chemical and biochemical engineering and scientific 
disciplines. For demonstration, the reaction-diffusion 
process occurring in a spherical porous immobilized 
biocatalyst has been solved successfully with the Michaelis- 
Menten form of kinetics. Simple analytical expressions of 
the concentrations and the effectiveness factors are derived 
by using the New Homotopy perturbation method. Our 
analytical results are compared with the numerical results 
and a satisfactory agreement is noted. 


[Parshotam et. al (1991)]. Another method based on a 
Taylor series linearization technique was given for a 
biofilm reactor with Monod kinetics and product 
inhibition [Bhaskar et. al (1991), Walas et. al (1995) and 
Silva et. al (2001)]. 

In this paper, a simple analytical expression is 
suggested to identify the substrate concentration 
distribution inside the immobilized biocatalyst 
particles and further to estimate their effectiveness 
factor. This offers an option for fast estimation of such 
a parameter in practice. 
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Introduction 

The rising interest in the application of biocatalysts 
entrapped in inert porous supports is mainly 
ascribable to their high volumetric productivities and 
the characteristics of being used repeatedly. The 
effectiveness factor of these immobilized biocatalyst 
particles is of capital importance in optimizing the 
processes because usually the diffusion of the 
substrates into the livelihoods and the products out of 
the supports have been affected significantly. Many 
works have been carried out to create mathematical 
models using the reaction-diffusion equations, which 
are commonly solved by numerical means (of iterative 
nature) as there are limited analytical solutions [Li et. 
al.(2004) and Villadsen et. al (1978)]. On that point are 
also iterative approximate solutions available in 
literature, which are obtained by comparison theorems 


Mathematical Formulation Of The Boundary 
Value Problem 


The immobilized spherical biocatalyst particle of the 
mass transfer is ordered by the Fisk's first law and the 
biochemical reaction follows the Michaelis -Menten 
kinetics, the reaction-diffusion process of substrate in 
the particle at steady state may be written every bit: 


d 2 C 2 dC r C 

1 _ HL 

dr 2 r dr D(K m +C ) 


(1) 


To reduce the eqn.(l) into the normalized form, we use 
the following dimensionless variables: 






( 2 ) 


Now, the eqn. (1) become 


±JL + l*L- 9q ? 

dx x dx 


}+py. 


= o 


(3) 


The corresponding boundary conditions of the eqn. (3) 
are 
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y = 1 at x = 1 

(5) 

The second boundary condition assumes 

the external 

convection resistance is neglected. By the definition of 
the effectiveness factor for a spherical immobilized 

biocatalyst, an equation of it is 


47rR 2 D(dC/ dr} r _ R 
(4/3) 7rR 3 (r m C s /(K m+ C s )) 

(6) 

B + \( dy\ 

J7 = s — r \ — 

3 <P \dx ) x= i 

(7) 


Analytical Expression Of The Normalized 
Surface Concentrations Using The New 
Homotopy Perturbation Method 

Linear and non-linear phenomena are of fundamental 
importance in various fields of science and 
engineering. Most models of real life problems are still 
very difficult to solve. Therefore, approximate 
analytical solutions such as a Homotopyperturbation 
method (HPM) [Gohri et. al (2007), Ozis et. al (2007), 
Li et. al (2006), Mousa et. al (2008), He (1999 and 
2003), Ariel (2010), Ananthaswamy et. al (2012 and 
2013) and Shanthi et. al (2013)] were introduced. This 
method is the most effective and convenient one for 
both linear and non-linear equations. Perturbation 
method is based on assuming a small parameter. The 
majority of non-linear problems, especially those 
having strong non-linearity, have no small parameters 
at all and the approximate solutions obtained by the 
perturbation methods, in most cases, are valid only for 
small values of the small parameter. Generally, the 
perturbation solutions are uniformly valid as long as a 
scientific system parameter is small. However, we 
cannot rely fully on the approximations, because there 
is no criterion on which the small parameter should 
exist. Thus, it is essential to check the validity of the 
approximations numerically and/or experimentally. 
To overcome these difficulties, HPM have been 
proposed recently. 

Recently, many authors have applied the Homotopy 
perturbation method (HPM) to solve the nonlinear 
boundary value problem in physics and engineering 
sciences [Gohri et. al (2007), Ozis et. al (2007), Li et. al 
(2006), Mousa et. al (2008)]. Recently this method is 
also used to solve some of the non-linear problem in 
physical sciences [Mousa et. al (2008), He (1999 and 
2003)]. This method is a combination of Homotopy in 


topology and classic perturbation techniques. Ji-Huan 
He used to solve the Lighthill equation [Mousa et. al 
(2008)], the Diffusion equation [He (1999)] and the 
Blasius equation [He (2003)]. The HPM is unique in its 
applicability, accuracy and efficiency. The HPM and 
New HPM [Gohri et. al (2007), Ozis et. al (2007), Li et. 
al (2006), Mousa et. al (2008), He (1999 and 2003), 
Ariel (2010), Ananthaswamy et. al (2012 and 2013) and 
Shanthi et. al (2013)] uses the imbedding parameter 
up as a small parameter, and only a few iterations are 
needed to search for an asymptotic solution. The 
approximate analytical expression of concentration 
(see Appendix B) of the icons. (3)-(5) using the New 
Homotopy perturbation method [Ananthaswamy et. 
al (2012 and 2013), Shanthi et. al (2013)] is given by 

sinh(£x) /ON 

y = (o) 

x sinh(£) 


The corresponding effectiveness factor rj can be 
calculated by using the eqn. (7) is as follows: 


/7 = 


Qg + 1) 

3 (p 1 sinh(&) 


(k cosh (k) + sinh (k)) 


(9) 


where k is defined by 


k - 


3 (p 

V(i+/?) 


(10) 


Numerical Simulation 

The non-linear reaction-diffusion equations (3 )-(5) are 
also solved by numerically. We have used the function 
pdex4 in Scilab/Matlab numerical software to solve 
numerically, the initial-boundary value problems for 
parabolic-elliptic partial differential equations. This 
numerical solution is compared with our analytical 
results (Figs. (l)-(2)). 

Results And Discussions 

Figured shows the dimensionless substrate 
concentration y{x) versus the dimensionless radial 
distance x - From these Figs, it is clear that when the 
dimensionless Michaeliss-Menten constant /? 
increases, the corresponding substrate concentration 
y(x) also increases for various values of the Thiele 

modulus cp 1 . Figure.2 also shows the dimensionless 
substrate concentration y{x) versus the dimensionless 
radial distance x - From these Figs, we infer that when 
the Thiele modulus cp 1 increases, the corresponding 
dimensionless substrate concentration y(x) decreases 
for various values of the dimensionless Michaeliss- 
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Menten constant /? . 

Figure. 3 shows the dimensionless effectivenss factors 
77 versus the dimensionless Michaeliss-Menten 
constant /? . From these Figs, it is observed that when 
the Thiele modulus cp 2 increases, the corresponding 
effectiveness factors 77 also increases. Figure. 4 , shows 
the dimensionless effectivenss factors 77 versus the 

(a) 
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(d) 



FIGURE 1 DIMENSIONLESS SUBSTRATE CONCENTRATION y 
VERSUS DIMENSION RADIAL DISTANCE x COMPUTED BY 
USING THE EQNS. (8) AND (10) FOR ALL VALUES OF THE 
DIMENSIONLESS MICHAELISS-MENTEN CONSTANT p, AND 

FOR VARIOUS VALUES OF THE THIELE MODULUS <p 2 , 
WHEN (A) <p 2 = 0.01,(B) (p 2 = 0.05, (C) (p 2 = 1 AND (D) (p 2 = 2. 

(a) 
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(c) 



FIGURE 2 DIMENSIONLESS SUBSTRATE CONCENTRATION y 
VARIOUS DIMENSION RADIAL DISTANCE * COMPUTED BY 
USING THE EQNS. (8) AND (10) FOR ALL VALUES OF THE OF 
THE THIELE MODULUS cp 2 , AND FOR VARIOUS VALUES 
DIMENSIONLESS MICHAELISS-MENTEN CONSTANT fl, WHEN 
(A) >0=0.01, (B) p= 0.1 AND (C) >0 = 0.2. 


(a) 



(b) 



FIGURE 3 DIMENSIONLESS MICHAELISS-MENTEN CONSTANT 
P VERSUS THE EFFECTIVENESS FACTOR r/ FOR VARIOUS 

VALUES OF THE THIELE MODULUS cp 2 , WHEN 
(A) (p 2 = 2,4,6,8,10,12 and 14 (B) p 2 = 0.01, 0.1 and0.2. 
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(a) 




FIGURE 4 THIELE MODULUS cp 2 VERSUS THE EFFECTIVENESS 
FACTOR ;; FOR VARIOUS VALUES OF THE MICHAELISS- 
MENTEN CONSTANT p , WHEN (A) p = 2, 4, 6, 8, 10, 12 and 14 , (B) 
>0 = 0.01,0.1 and 0.2. 

dimensionless the Thiele modulus cp 2 . From these 
Figs, we note that when the dimensionless Michaeliss- 
Menten constant p increases, the corresponding 
effectivenss factors 77 also increases. 

Conclusions 

This research article reports a mathematical treatment 
for immobilized spherical biocatalyst of the mass 
transfer. The novelty of this manuscript is an 
application of approximate method of the second 
order non-linear partial differential equations. The 
approximate analytical expressions for the steady state 
substrate concentration profiles for all values of the 
dimensionless parameters p and cp 2 are obtained 
using the New Homotopy perturbation method. 
Furthermore, an analytical expression corresponding 
to the steady state effectiveness factor is also 
presented. The effectiveness factors of the biocatalyst 
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particles may be estimated using the HPM for some 
values of cp 2 such as oxygen transport into natural 
my celia pellets in submerged cultures. A satisfactory 
agreement with the numerical results is noted. This 
analytical expression presented in this manuscript can 
be easily extended too many non-linear reaction- 
diffusion equations in biosensors. 


first use the embedding parameter as a small 
parameter and assume that the solutions of the eqns. 
(A.4) and (A. 5) can be written as a power series in p : 

v = v 0 + pv x + p 2 v 2 + ... (A. 8) 

Setting p = 1 results in the approximate solution of 
the eqn. (A.l): 
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Appendix A: Basic Concept Of The 
Homotopy Perturbation Method 

To explain this method, let us consider the following 
function: 

D 0 (u)-f(r) = 0, reQ (A.l) 

with the boundary conditions of 

B a (u, — ) = 0, rer (A.l) 

on 

where D 0 is a general differential operator, B 0 is a 
boundary operator, / (r) is a known analytical 
function and T is the boundary of the domain Cl . In 
general, the operator D 0 can be divided into a linear 
part L and a non-linear part N . The eqn. (A.l) can 
therefore be written as 

L(u) + N(u)-f(r) = 0 (A.3) 

By the Homotopy technique, we construct a 
Homotopy v(r, p ) : Cl x [0, 1] — » 91 that satisfies 

H(y , p) = (1 - p)[L(y) - L(u 0 )] + p[D 0 (v) - /(r)] = 0 (A.4) 
H(v , p) = L(y) — L(u 0 ) + P L(u 0 ) + p[N(y) - f(r)\ = 0 (A.5) 

where p e [0, 1] is an embedding parameter, and u 0 is 
an initial approximation of the eqn. (A.l) that satisfies 
the boundary conditions. From the eqns. (A.4) and 
(A.5), we have 

H(v, 0) = L(y) - L(u 0 ) = 0 (A.6) 

H (v, 1) = D 0 (v) - f(r ) = 0 (A. 7) 

When p= 0, the eqns. (A.4) and (A.5) become linear 
equations. When p =1, they become non-linear 
equations. The process of changing p from zero to 
unity is that of L(y) - L(u 0 ) = 0toD o (y)-f(r) = 0. We 


u - lim v = v 0 + Vj + v 2 + ... 
This is the basic idea of the HPM. 


Appendix B: Analytical Solution The Eqns. 
(3)-(5) Using The New Homotopy 
Perturbation Method 

In this Appendix, we indicate how the eqn. (8) in this 
paper is derived. To find the solution of eqns.(3) - (5) 
we construct the new Homotopy as follows [shanthi 
et. al (2013)] 

( !_ jP) fe + 2 tevr 9^] J +/ ,fe + 2 r^_r = 0 

dx pydx) ll + pj dx p\dx) ll + /?y) 


The analytical solution of the eqn.(B.l) is 


T^To+PTi+P y 2 +- 


Substituting the eqn.(B.2) into an eqn. (B.l) we get 


d [yo+pyi+p yi+- 


d (yo+pyi+p ^ 2 +•• 


, 2 dly 0 + py { + p 7 2 + - 2 d(yo + M + P 72 + - 

1 -p) + +p + =0 

x dx x dx 


V l 2 
(yo + py\ + p 72 +•••) 


9 ^ (y 0 + m + P 72 + -) 
1 + A7o+/ ? 7i+/72+-) 


Comparing the coefficients of like powers of p an 
eqn.(B.3) we get 

< b - 4 > 

dx p\ dx ) ^ 1 + p ) 

The initial approximations are as follows: 

y'o (0) = 0 and y 0 (i) = l (B.5) 

7(0) = 0 and >-,(1) = for i = 1,2,3 (B.6) 

Solving the eqn. (B.4) and using the boundary 
conditions the eqns.(B.5)-(B.6), we obtain the following 
results: 

y = (B.8) 

x sinh(A:) 
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where k is defined by the text eqn. (10) 

According to the HPM, we can conclude that 

y = \imy(x) = y 0 (B.9) 

p-+l 

After putting the eqn.(B.8) into an eqn.(B.9), we obtain 
the solution in the text eqn. (8). 

Appendix C: Matlab/Scilab Program To Find 
The Solution For The Non-Linear Eqns. 
( 3 )-( 5 ) 

function 
m = 2; 

x = linspace(0,l); 
t = linspace(0,100000); 

sol = pdepe(m,@pdex4pde,@pdex4ic,@pdex4bc,x,t); 

ul = sol(:,:,l); 

figure 

plot(x,ul(end,:)) 
titleful (x,t)') 
xlabelfDistance x') 
ylabel('ul(x,l)') 

function [c,f,s] = pdex4pde(x,t,u,DuDx) 
c =1; 

f=l * DuDx; 
phi=sqr t (0 . 05 ) ;be ta=l 4 ; 

F=-9^(phi) A 2 >f u(l)/(l+(beta >f u(l))); 
s =F; 

function uO = pdex4ic(x); 
uO = 0; 

function [pl,ql,pr,qr] = pdex4bc(xl,ul,xr,ur,t) 
pi =0; 

q i = i; 

pr = ur(l)-l; 
qr = 0; 


Appendix D: Nomenclature 


Symbol 

Meaning 

y 

Dimensionless substrate concentration 

p 

Dimensionless Michaelis Menten constant 

V 

Effectiveness factor 

X 

Dimensionless radial distance 

V 

Thiele modulus 

c 

Substrate concentration 

C s 

Substrate concentration on the surface of the 
particle 

D 

Effective diffusivity 

K m 

Michaelis-Menten constant 

r 

Radial distance 

r m 

Maximum reaction rate' 

R 

Particle radius 
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